	subroutine shelf_qs
	include 'cst.inc'
c
c   calculates the time-, space-, & depth-averaged sediment
c   transport rates. k1 and beta control the
c   net onshore sediment flux due to the waves. covers the whole
c   shelf and shoreface but not the surf zone.
c	 
c	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c	3 ,save1(200),save2(200)
c	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c
c  !!!!!!!!!!!!!!!!!!!!!!!!awn5/27!!!!!!!!!!!!!!!!!!!!!
c	do 774 jxx=1,25
c	j=jxx
c774  write(44,*)'j =',j,'cstar =',cstar(25,j),'h=',h(25,j)

c  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c

	do i=1,imax
	   j=jshore(i)-1 !Exclude the shoreline cells
c	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!
c	if(j.lt.20)write(*,*)'shelfqs i =',i,'     j=',j
c	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         slope(i) = (h(i,j+1)-h(i,j))/deltay
      enddo   
c
c
	do i=1,imax
	   do j=1,jmax  !30 needs to be passed thorugh as jmax
	      qyup(i,j) = 0.0 !Upper face facing north
	   enddo   
      enddo
c
	ilim=imax-1
	do i=1,ilim
	   jlim=jshore(i)-1		  !was originally jshore(i)-2
	   do j=1,jlim
c
c   determine the values for the upper face of each cell
c   
	      dlup = (dl(i,j)+dl(i,j+1))/2. !Diffusion coefficient of long-shore direction
	      dcup = (dc(i,j)+dc(i,j+1))/2. !                         cross-shore
	      aup = (a(i,j)+a(i,j+1))/2.    !used only for correction
	      awup = (aw(i,j)+aw(i,j+1))/2. !used only for correction
	      hup = (h(i,j)+h(i,j+1))/2.
c
c    determine the y-diffusion coef. from the long-, and x-shore values
c
	      dybup = dcup !*(cos(aup))+dlup*(sin(aup))			   !change later
c
c    determine the diffusive and wave-driven transport for the up face
c    and total these
c??????????????????????????????????The angle correction terms are very important 
            !sediment flux of diffusion along the y-direction
	      qy1up = -k1*dybup*((cstar(i,j+1)-cstar(i,j))/deltay)
	      qy2up = k1*wo*exp(beta*hup) !*(cos(awup))			  !change later
	      qyup(i,j) = qy1up + qy2up
         enddo
      enddo
c
	do i=1,imax
	   do j=1,jmax  !30 needs to be passed thorugh as jmax
	      qxrt(i,j) = 0.0
         enddo
      enddo   

	ilim=imax-1
	do i=1,ilim
	   jlim=jshore(i)		  
	   do j=1,jlim
c
c    determine the values for the right face of each cell
c
	      dlrt = (dl(i,j)+dl(i+1,j))/2.
	      dcrt = (dc(i,j)+dc(i+1,j))/2.
	      art = (a(i,j)+a(i+1,j))/2.
	      awrt = (aw(i,j)+aw(i+1,j))/2.
	      hrt = (h(i,j)+h(i+1,j))/2.
c
c    determine the x-diffusion coef. from the long-, and x-shore values
c
	      dxbrt = dlrt*(cos(art))+dcrt*(sin(art))
c
c    determine the diffusive and wave-driven transport for the rt face
c    and total these
c
	      qx1rt = k2*dxbrt*((cstar(i,j)-cstar(i+1,j))/deltax)
	      qx2rt = k1*wo*exp(beta*hrt) !*(sin(awrt))
	      qxrt(i,j) = qx1rt + qx2rt
	   enddo   
	enddo   
c***************************All the lines below are not used	
c
c	now modify qxrt values for interaction with jshore cells
c     This is to check whether there is flow in the x-direction in shoreline cells
c
	do i=1,imax-1
	   js1 = jshore(i)    !index of one cell
	   js2 = jshore(i+1)  !index of the next cell
	   if(js1.eq.js2) qxrt(i,js1) = 0.0	!If the two cells are in the same horizon, there is no flow between them
	   if(js1-js2.eq.1) qxrt(i,js1)=0.0   !JS1 is northwest to JS2, you cannot put anything to the land
	   !This if condition is never met
	   if(js1-js2.eq.2) then
	     qxrt(i,js1) = 0.0
	     qxrt(i,js1-1) = 0.0
	   endif
	   if(js1-js2.eq.-1) then
	   endif
	   if(js1-js2.eq.-2) then
	   endif
	   !End here
	   if(abs(js1-js2).gt.2) then
	      write(*,*)'jshore variable differs by more than 1 cell'
	      write(*,*)'between row ',i,' and ',i+1
	      write(*,*)'jshore(',i,')= ',jshore(i)
	      write(*,*)'jshore(',i+1,')= ',jshore(i+1)
	      stop
	   endif
	enddo
c
c	ordinarily the sed. flux on the landward (-up)
c	face of the surf zone cell (jshore(i)) is zero.
c	this whole set of cells can migrate bodily as
c	they empty or fill.
c	the transition cells represent the upper shoreface
c	which is taken to have a time=averaged constant profile
c	characterized by a slope (scr).
c
c
c
c   there are special conditions in the surf zone and transition cells.
c   consider first the cells in the transition zone, just seaward
c   of the surf zone.  we take the cross shore flux (qy) across its
c   seaward face in the same way that we do all of the other shelf
c   cells.  the flux across its face shared with the surf zone cell is
c   taken to be proportional to the difference of the local bottom slope
c   and a critial value of this slope (scr).  if this slope is larger than a
c   given value (scr) then qy across this face is offshore (negative). 
c    this flux is removed directly from the surf zone.
c
c
c   set left and right constant boundary conditions
c
cwr	jlim=jshore(1)
cwr	do 50 j=1,jlim
cwr	qxrt(1,j) = qxrt(2,j)
cwr	qyup(1,j) = qyup(2,j)
cwr   50	 continue
c
cwr	jlim=jshore(imax)
cwr	do 55 j=1,jlim
cwr	  qxrt(imax,j) = qxrt(imax-1,j)
cwr	  qyup(imax,j) = qyup(imax-1,j)
cwr   55	 continue
c
c   set grid lower boundary constant flux condition
c
cwr	do 60 i=1,imax
cwr	qxrt(i,1)=qxrt(i,2)
cwr   60	qyup(i,1)=qyup(i,2)


cccccccccccccccccccccccccccccc
c	do 789 i=1,imax
c	i=5
c	do 789 j=1,jshore(i)+2
c	write(88,*)i,j,qyup(i,j)
c	qxrt(i,j) = 0.
c789	continue

	return
	end